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We study the consequences of noise and dissipation for parametric resonance during preheating. 
The effective equations of motion for the inflaton and the radiation field are obtained and shown 
to present self-consistent noise and dissipation terms. The equations exhibit the usual parametric 
resonance phenomenon, allowing for exponential amplification of the radiation modes inside the 
instability bands. By focusing on the dimension of the border of those bands we explicitly show 
that they are fractal, indicating the strong dependence of the outcome in the initial conditions. The 
simultaneous effect of noise and dissipation to the fractality of the borders are then examined. 
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I. INTRODUCTION 



In standard inflationary scenarios of the early universe, the transition from the inflationary regime itself to the 
radiation dominated phase (or standard big bang regime) is accomplished through a short period of reheating after 
inflation, when the inflaton field, by rapidly oscillating at the bottom of its potential in a time scale shorter than 
J> the Hubble time, releases all of its energy in the form of light particles that it may be coupled to. Recent studies 
T-H , on this reheating phase have shown the possibility of emergence of new and interesting phenomena, like explosive 
' particle production, nonequilibrium symmetry restoration, among others effects . These effects are directly related 
" to the possibility of the oscillating inflaton fleld to drive the fluctuations of the coupled flelds to regimes of parametric 
ampliflcation, or parametric resonance, where particles are efficiently produced. 

The parametric resonance, or preheating phase, has then been shown to be a fundamental effect that must take 
f-"^ . place before the actual reheating of the universe through the thermalization of the produced particles ||^ . Preheating 
also can be easily obtained in very different models by means of a linear or a quadratic coupling between the inflaton 
and the radiation fields. In all those cases the inflaton, oscillating at the bottom of its potential (after the inflationary 
I ■ phase), behaves as a periodic mass for the other scalar fields that are coupled to it. As a consequence, the equations 
Qh, for the modes k of the fluctuation fields present exponential solutions 
D . 
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Xfe(0 ~exp[±/i-t] , (1.1) 



where /i, the Floquet exponent, may be either real or complex, depending on the parameters adopted. One is thus able 
^ I to define the so called instability bands in the parameter space, within which the Floquet exponent is real, indicating 
■ - ' exponentially (de)increasing solutions. The regions of the parameter space where /i is imaginary, on the other hand, 
correspond to oscillatory solutions and are thus called stability bands. 

One important aspect concerning the band structure and exponential growth of fluctuations is the possible chaotic 
behavior associated to the very existence of the bands, both at their borders and inside them. Our concern in this 
paper is thus to determine the dimension of the borders between two such bands. We will also be interested in 
calculating the dimension of the border between two different values of the Floquet exponent inside a given instability 
band. As we will show, all of them are fractal, unambiguously characterizing them as chaotic ones. As such, they 
prevent one from telling if a given mode, assumed close enough to a border, will be amplifled or not. Thus, the 
outcome of the preheating phase will strongly depend on the state of the universe right after the inflationary period. 

The fractality, or chaoticity, of the borders separating instability bands is, in a sense, not a completely unexpected 
result. For example, in a recent work ||^ the authors argued for the close relation between the presence of chaos and 
parametric resonance. In fact, the equation of motion for the fluctuation fields coupled to the inflaton field, viewed 
as a classical equation of motion, resembles that of an harmonic oscillator driven by a periodic external force, which 
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is well known to exhibit chaotic behavior under appropriate conditions. A few other papers have also studied chaos 
in the context of preheating, as by the authors of Ref. |Q, although they have studied chaos after the preheating 
period. In general, the dynamics of interacting fields may show the possibility of chaotic dynamics, as has been shown 
recently in Refs. g ^. 

One important aspect associated to the microscopic physics of interacting fields raised in Refs. |j was the 
importance of dissipation in the dynamics of coupled system of fields and the influence of dissipation on the degree 
of chaoticity of the field dynamics. Recall that typically the effect of dissipation is to damp the fluctuations on the 
system and consequently tends to suppress possible chaotic motions and makes field trajectories in phase space to 
tend faster to the system asymptotic states. This is indeed the case for the dynamics of background (or classical) field 
configurations. However, the fluctuation (or quantum) modes of the fields that experience parametric resonance are 
not only subject to dissipation but also to stochastic noise, being both related by a generalized fluctuation-dissipation 
relation 0. 

Previous works have also studied the influence of noise during preheating. For example the authors of Refs. ||] 
have studied the consequences to preheating when a somewhat ad hoc noise is added as a mass term to the secondary 
field's equation of motion. Here, in contrast, we will take the equation for the inflaton as a typical ensemble averaged 
equation of motion (as we should expect for a background field configuration) and noise is only included in the effective 
equation of motion for the fluctuation field coupled to the inflaton. This equation is self-consistently derived for both 
dissipation and noise through the usual method of Schwinger's closed time path formalism of real time (for previous 
references, see for instance Refs. [Q, 0). The simultaneous effect of noise and dissipation for the reheating 
have been also considered recently by the authors of Ref. |]l2j where they also derived self-consistently the effective 
equations of motions for both the inflaton and the radiation fields. Their main interest was, however, in the final 
stages of reheating, not including therefore the consequences of both noise and dissipation for the preheating phase. 
This is our main concern in this paper, which we will specialize in the simultaneous effect of noise and dissipation on 
the fractality of the borders of the instability bands and how these microscopic physics associated to the dynamics of 
the fields may alter the parametric resonance phenomenon during preheating. 

The chaotic behavior associated to the parametric resonance stage is then studied by means of the determination 
of the fractal dimensions of the borders of the instability bands. As argued before in Refs. |5[ ^ (see also |l^), the 
fractal dimension gives a topological measure of chaos for different space-time settings and it is a quantity invariant 
under coordinate transformations, yielding then an unambiguous signal for chaos in cosmology and general relativity 
problems in general. The dimension of the border provides thus a robust characterization for chaos as opposed to 
Lyapunov exponents which may change sign by a simple change of observer p^ . 

This paper is organized as follows: in Sec. H we consider a simple model leading to preheating consisting of an 
inflaton field (j> of chaotic^ inflation with quadratic potential coupled to a set of radiation scalar fields Xj that cj) can 
decay to. Additionally we allow the radiation fields to decay to fermion fields ipk (m this model the final radiation 
energy density is composed of fermion matter fields). The effective equations of motion for both the inflaton (f> field 
and the radiation Xj fields are then explicitly derived. In Sec. HI we present our numerical studies concerning the 
parametric resonance in our model and the computation of the fractal dimension of the border of the instability bands. 
Finally in Sec. IV we give our concluding remarks. 

II. DERIVING THE EFFECTIVE EQUATION OF MOTIONS 

We will study here the simplest model of an inflaton field (p coupled to a set of scalar fields Xji 3 — 1, . ■ . , A^x' 
through a trilinear coupling. Although quadratic couplings between the fields are more common in the study of 
resonance during preheating, models with trilinear couplings exhibit the same phenomenon as well and had been used 
before as a toy model in different analysis ||l5| . We choose this simplest model here for convenience since it will allow 
a more straightforward analysis of the problem without having to compute higher order terms, as we are going to 
see below. Additionally, we couple to the radiation scalar fields Xj ^ set of fermion fields ^k, k = 1, . . . ,N^, thus 
allowing Xj to decay into the fermion fields. These kind of models as this one we will study here, allowing cascade 
decay sequences as </) — > 2xj and Xj ~^ tpk + 'tpk have shown to exhibit interesting dissipation properties and would 
then be particularly interesting to study them in the context of reheating and parametric resonance after inflation. 
The model we study here is then described by the following Lagrangian density 



Here, "chaotic inflation" only means we use a single-well potential for the inflaton, and has no relation whatsoever to the actual chaotic 
behavior we describe later on. 
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(2.1) 



A. The effective equation of motion for the inflaton field 



We can easily obtain the exact equation of motion for the inflaton field fi'om the model Eq. {2A). Decomposing 
the field into a classical background component ip (the classical inflaton field) and a quantum fiuctuation part as 
(j) = ip + Scf), with {4>) = if and (Scf)) = 0, then the equation of motion for ip readily follows by imposing {dcj)) — at all 
orders in perturbation theory (the tadpole method). For a homogeneous inflaton field this gives 
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^(t) + ™|^(t) + g|(x,^) = o 



(2.2) 



By perturbatively expanding (xj) in the ip amplitude, we obtain [ p!6| 



(2.3) 



where Gj+(q, t — t') is the real time causal propagator (in momentum-space) for the x held. The expression for 
G'J^^(q, t — t') can be found for example in Refs. [0, |l6[ p7| . 

Eq. ( |2.3| ) can also be directly obtained from the standard rules of the real time formalism. In the absence of 
fermions, we can easily functionally integrate the Xj held exactly, since Eq. (2.1) has only quadratic terms, to obtain 



r[^] = S[4>] + -trln[n + + gjcj,] , 



where 



-Trlnp + m^^ + g'^<P] = -iln / I?Xj exp j --Xj [□ + "^x, + 9j<P\X3 



(2.4) 



(2.5) 



and the field is integrated over a time path c that goes from —00 to +00 and then back to —00. This is the Schwinger's 
closed time path formalism. 



Eq. (2.4) is a nonlocal equation in <f). By perturbatively expanding (2.4) in the inflaton amplitude (j) (or in the 
coupling constant g, for g ^ 1) we obtain (see for example Ref. H) 



-Trln [u + m\+g'^c, 



(-1)'"+! 



>^ E ^^G™ {gS 



{g^Y' Tr j d^x, . . . d^x^Gl^''^ (xi - X2) ['/'(x2)]i^,„, Gl-''-{x2 ~ x^) . . 



[<^(a^™)]u-i,n™ G;;"''"(a;„ - xi) [0(xi)];^ „^^^ , 



where G"''(a; — a;') is the real time x-held propagator on the contour c, given by {l,n = +, — ) ||Tq, 



(2.6) 



G++ix-x')^i{T+xix)xix')) 

G--{x-x')=i{T^x{x)x{^')) 
Gj"(a; - x ) = i{x{x')x{x)) 
G~+{x-x')=i{x{x)x{x')), 



(2.7) 
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where T+ and T_ indicate chronological and anti-chronological ordering, respectively. is the usual physical 

(causal) propagator. The other three propagators come as a consequence of the time contour and are considered as 
auxiliary (unphysical) propagators The explicit expressions for G^\x — x') in terms of its momentum space 

Fourier transforms are given by 20| 



\G- + {ci,t-t') G--{q,t-t') 



(2.8) 



where 



G++ (q, t~t') = G> (q, t - t')0{t - t') + G< (q, t~t')eit'-t) 
G^- (q, t-t') = G> (q, t - t')e{t' -t) + G< (q, t - t')e{t - t') 
G+-{q,t-t') = G<{q,t-t') 
G-+{q,t-t') = G>{q,t~t') 



(2.9) 



and (we formally here consider the full propagator expression for the x held. This will be convenient later, and we use 
the approximation that the spectral function for x has a Breit-Wigner form, with width , to write the expressions 
for the propagator as |l^ ) 



G>{q,t-t') 



1 

2^ 



{exp[-i(tJx - ■i^x)it - - t') + exp[-i{uj^ + iT^){t - t')]9{t' - t)} 



G<iq,t-t')=G>iq,t' -t) 



(2.10) 



For the analogous expressions at finite temperature, see for example Ref. pq |. 

In terms of the field variables 0+ and 4>- (in the forward and backward branches of the time contour, respectively) 
and the real time propagator Eq. (p^), the effective action to order becomes 



(5^0-) 



2 "^0^2 



2^ t^-r y 2 

+ {x)G++{q,0)-cf,-{x)G-~{q,0)] 



+ i 



^ / d^xd^-' 



(2^)3 Lv+^-^^^x y^^-' -^-v-y-x 
d^k , r d^q 



4 J 7 (|?^^^P[*^-(^-^')] / ^[Mx)G++{k + q,t-t')G++iq,t~t')M^') 

- M^)G+- (k + q, i - t')G+- (q, t ~ t')(t>^ {x') - {x)G-+{k + q,t~ t')G-+{q, t - t')cp+{x') 

+ ^^{x)G--{\i + q,t~t')G--{q,t-t')^^{x')\ . (2.11) 

At this point it is more convenient to introduce two new variables 0c and <j)/^, given in terms of 0+ and by 



(2.12) 



Rewriting the effective action ( 2.11| ) in terms of ipc and i/ia, we obtain 



/g^ f f d^q 1 

d'^x(t)A{x)[~n ~ ml](t)c{x) - Y J d'^xcjjAix) J J^^^^ 

+ 1^1 dWa;VA(x)0A(x') J ^exp[*k.(x-x')] J ^Rc[G++ {k + q,t - t')G++ {q,t - t')] 



g''' J d'''xd'''x'(l)A{x)(pc{x') J ^^^^ exp[ik.(x - x')] 
/ (SyJ^^^t^^^^^ + ^ " i')G++(q, t - t')]e{t - t') 



(2.13) 



5 



The imaginary term in Eq. ( ^.13 ) can be associated as coming from a functional integration over a Gaussian noise 
(stochastic) field 



D^^P[^^]exp J d'^x(l)A{x)^4,{x) 

= exp|-y d^xd''x'^<PAix)Re[G++G++]^^^,Mx')^ , 
where P[^0] is the probability distribution for and it is given by 

P[e^] = Ar-iexp|-i J d4:rdVe^(x)(^^Re[Gj+Gj+]^ ^,) ^^x')^ 
In Eqs. ( 2.14 ) and ( ^.15| ) we have used the short notation for the product of propagators: 



(2.14) 



(2.15) 



_g_exp[zk.(x-x')] / ^[Gj+(k + q,t-t')Gr(q,t-0]- 



Using Eq. ( 2.14 ) then we can write Eq. ( 2.13 ) as 



ri^c^A] = -«ln J D^^P[^^]exp{iSes[(l)c:4'A,i,p]} , 



where 



S'cff[<^c,0A,^0] J d'^x(l)Aix)[-D - ■ml](f>c{x) -Y^^ — J d^xcjiAix) J — 



d^q 1 



)3 2a;, 



(2.16) 



(2.17) 



- Y.f^j d^xd^x'^A{x)Mx')l^ K^GJ^].,.' ^(^ - + / rf*^'/'A(x)C0(:E) . (2.18) 

The second term in the above equation is divergent. It can be easily removed by adding the appropriate renormalization 
counter-term in the classical potential. In the following we assume a renormalized action and just drop that divergent 
term from the equations. 

The equation of motion for 0c is defined by (see for instance Ref. 0) 



SScg[4'A,4>c,£,4>] 



from which we obtain 



(2.19) 



{D + ml)Mx)+Y.9' / d^x'Mx')lm [G++G++] 9{t - t') ^ ^^{x) . 



Note, from Eqs. (2.14) and (2.15), that the noise field satisfies 



(2.20) 



4 

(e4x)c4x')) = EfRe 



(2.21) 
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Taking the (ensemble) average of Eq. ( ^.2C ) and associating {4>c) to the classical inflaton field configuration, we 
obtain (for a homogeneous field, ip = (p(t)) 



This equation is exactly the same as the one obtained from the tadpole method, Eq. (2.3) 
Using Eqs. (2.9) and (2.10), we obtain that 



(2.22) 



Im[Gj+(q,t-t')] 



t>t' 



exp(-2r^Jt-<'|) . 

— sm{2uj^^\t — t \) 



Substituting this into Eq. (2.22), we obtain 



(2.23) 



d^q exp(-2rxJ<-i'|) . 

^2 sn.^.u.^^. 



sin(2w^.|i - t'\) ^ . 



Xi 



By integrating by parts in the time the last term in the above, we obtain 



(2.24) 



J — oc J 



d^q exp{-2r^^ \t - t'\) [uj^^ cos(2wx, 1^ - ^'D + ^x, sin(2cOx^. \t - t'\)] 



(27r)3 



where 



(2.25) 



ml = ml 



1 



,.1 ^ (2-)^8-x.(r^.+<) ■ 

Using now that the Xj fields are allowed to decay into the ip^ fields with a decay rate T^. (q) given by |lf 



(2.26) 



^xAq) 



3/2 



1-4- 



X] 



Xk 



(2.27) 



where a^. = J2k=i K'^lji'^ - ^'^U'^lJ^^'^ / i^^)- Substituting this into Eq. ( |2.25D and approximating the dissipa- 



tive kernel by a Markovian one (this is a valid approximation provided we have a sufficiently large number of decay 
channels available, see [|l6|) we find the local equation of motion for tp 



Nx 



dh 



^xM) 



(2.28) 



Substituting Eq. ( ^.27| ) in the above equation, the momentum integral in Eq. (2.28) can be performed and the final 
result for the effective EOM for the homogeneous field f. 



(p{t) + ml(p{t) + Tj^ytpit) = 



(2.29) 



where 



4 9 



,jtl287r(m4^.+Oi(2y<T<+2m2^.)^ ' 
This equation has a simple damped solution for Lp(t) (overdamped for 77^ > 2?7i0, underdamped for 77^ < 2771^) 



(2.30) 
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B. The Equation of Motion for the Fluctuations Xj 



We can derive the effective dissipative and stochastic equation for the fluctu ation fields x in a similar way. The 
perturbative effective action for x for the model (2.1) is ^ (in analogy with Eq. (2.13)) 



+ z|-y d^xd^x'xA{x)xA{x')Re[G++G++U,, 
- 2/ [ d^xd^x'xA{x)Xc{x')lm[G++G++l,,,e{t-t') . 



(2.31) 



Again, we can interpret the imaginary term in the above equation as coming from a functional integration over a 
noise field 



J ^Cx-PKx] exp |i J d'^xxA{x)^xi^ 



exp 



d'^xd'^'x' ^XA{x)Ke 



ri++r'++ 



XAix') > , 



where probability distribution for and it is given by 



(2.32) 



= exp <^ -- / d^d'^x'^xix) /Re 



X 



x,x' 



ixi^') 



(2.33) 



Rewriting the action for x in terms of the noise field then the equation of motion which we obtain (analogous to 
Eq. (OO)) is 



{U + ml+ 9^ip)Xc{x) + 2/ / d'^x'xc{x')lm 



^X 



(2.34) 



In order to get an approximate analytical expression for the above equation of motion, we will assume only contri- 
butions to the effective action with zero external momentum (the linear response approximation) or nearly spatially 
homogeneous fields. This way we can handle the spatial nonlocality in Eq. (2.34) and obtain: 



{n + ml+ g^(p)xc{x) - 2/ / dt'xd^, t') 



f d^q exp[-(rx + r^)|<-i'| 



Integrating by parts in the time the last term in the Ihs of the above equation, we obtain 



sm[{uj^+uj^)\t -t'W • (2-35) 



{U+ml+ g\)xc{x) + 2.9* / dt'xd^, 



exp[-(rx + T^)\t - t'\] {(r^ + r^) sin[{uj^ + uj^)\t - t'\] + {uj^ + uj^) cos[(cJx + ^^)\t - t'\]} 



= i^{x).{2.m) 



Once again, supposing valid a Markovian approximation for the dissipative kernel in the above equation, we obtain 
for the effective EOM for the fluctuation field Xj the expression 



^ In the following we neglect the corrections of the fermions to Xji given by {'ipkM>k)- This is fine in the Markovian limit, in which case 
the contributions from the fermions to dissipation and noise vanish as a consequence of having r^j_ = [tq| and the fermions only give 
a renormalization to the Xj fields mass. 
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(2.37) 



where 



(2.38) 



while the effective (unrenormaHzed) mass m^^. appearing in Eq. (2.37) and adding the contribution coming from the 
fermions, we find that 



- 2 2 

^X, ^ ^X, 



(27r)3 ^^^uj^ [(wx, + u^4,Y + (r^, + r^)2 



k=l 



(2.39) 



For the model (2J) we have F^^. as given by Eq. ( 2.27 ) and F^ is given by 



/ ^2 ■ 



^ 167rw0 



ml 



1/2 



(2.40) 



The noise field in Eq. ( ^.37[ ), from Eqs. (|2.32D and (|2.33| ), obeys 



(ex(2^)ex(2^'))=ff|Re[Gj+G++],,,, . 



(2.41) 



III. THE (IN)STABILITY BANDS AND THEIR BORDERS 



The equations (|2.29D and ( |2.37| ), obtained in the previous sections, are the basic equations of this paper. If 
one neglects their dissipation and inhomogeneous terms, the system will present the exact parametric resonance 
phenomenon, with instability bands in the parameter space inside which modes grow indefinitely, with a real (and 
positive) Floquet exponent. Here, we make no such approximation. We only restrict ourselves to a time period when 
the back-reaction can be safely neglected, i.e., in the beginning of the preheating phase. 

Our goal is twofold. First, we check the very existence of the band structure when one uses the set of equations 
mentioned above, obtained from a microscopic calculation, as opposed to inserting ad hoc terms in their classical 
counterparts. As we will shortly see, the parametric resonance is barely changed by this modification. Second, we 
investigate the dimension of the borders of the bands. It has been known p3| that a fractal border prevents one 
from determining the final attractor for which a given trajectory will tend to. Although we cannot formally define 
an attractor in our problem, we can split the trajectories in two different cases: the one which are amplified and the 
ones which are not. We will explain the exact procedure used later on this section. 

The instability bands are manifest in the momentum space k for the Eq. ( ^.37 ). A ctua lly, we will plot the bands 



in the space k x g , as has become usual in the literature. In momentum space Eq. (2.37) becomes 



X,(k, t)+ + + .g2<^(t) xj(k, t) + ?7;,4,(k, t) = (k, t) 



(3.1) 



which together with Eq. (2.29) form a coupled system of differential equations that we study numerically. The typical 



values of parameters that we use are 



15 



9j 
hk 

ipit = 0) 



g = 10-4^4(2 
h ^ 0.1 
Mpi/5 
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X,(k,< = 0) « 1 1 Mil (3.2) 

3/7 X lO^^Mpi 
m.0 1/7 X \Q-^Mpi 

where we have expressed all diniensionfull quantities with respect to the Planck mass Mpi. 

The chaotic behavior of the dynamical system of equations, as determined by the equations of motion for (p{t), 
Eq. (|^), and the one for the fluctuation field modes, Xj(k, t), Eq. (|3.l[), is then quantified by means of the 



determination of the fractal dimension of the boundary separating each of the instabilities bands. To define exactly 
what is to be called an unstable band we establish a set of cutoff values for the energy E of the field: Eq, \OEq and 
IOO-Eqj where Eq is the largest possible initial energy (in momentum space) in the range of initial conditions taken, 



2 



x/(k,fo) + fc'x'(k,to) + m2 v'(k,io)+g>(to)x?(k,io) • (3.3) 



Fig. |l| shows the band structure in the phase space and the set of first three instability bands we will restrict ourselves 
to. Each energy cutoff defines the border between two different ranges for the Floquet exponent. 

We next move to the determination of the fractal dimension. The fractal dimension is associated with the possible 
different exit modes under small changes of the initial conditions at t = and it will give a measure of the degree of 
chaos of our dynamical system. The exit modes we refer to above are the ones associated to trajectories that evolve 
to the stable band or to the unstable one, as determined by the given cutoff energy E,,. The method we employ to 
determine the fractal dimension is the box-counting method, which is a standard method for determining the fractal 
dimension of boundaries [ pT| . Its definition and the specific numerical implementation we use here have been described 
in details in Ref. [|||. 

The basic procedure is: given a set of initial conditions Xq at t = t^ that leads to a certain outcome for the trajectory 
in phase space, we study whether there will be a change of outcome for the trajectory due to a perturbation 5 (i.e., 
whether the perturbation will lead to a different attractor or not). Given a volume region in phase space around a 
boundary between different attractors and perturbing a large set of initial conditions inside that region, the fraction 
of uncertain trajectories, /(£), that result in a different outcome under a small perturbation can be shown to scale 
with the perturbation 5 as [Ell f{5) ~ 5*^, where e is called the uncertainty exponent. The box-counting dimension of 
the boundary in phase space separating different attractors, or fractal dimension /d, is given by ||2l|| fd = d — e, where 
d denotes the dimension of the phase space. For our system of equations of motion, d = A. For a fractal boundary 
fd>d—l, implying that e < 1, whereas for a non- fractal boundary, = d — 1, and e = 1. 

Following the method of box-counting for this problem, it is enough to consider variations only on the parameters fc^ 
and (7^, the wave-number squared and coupling constant squared, respectively. Perturbations around these variables 
are taken between 10~^ to 10"'', as we move from the first to the third band in Fig. 1 and where the initial conditions 
are determined. For each region separating the instability bands, set by the variations of k'^ and g^, we then consider 
a large number of random points (a total of 10.000 random points for each run was used) enough to produce a reliable 
statistics with errors around the percent (see Ref. |Q for a description of estimates of the statistical error) . All initial 
conditions are then numerically evolved by using an eighth-order Runge-Kutta integration method and the fractal 
dimension is obtained by statistically studying the outcome of each initial condition for each run of the large set of 
points. 

The fractal dimension of each of the borders shown in Fig. 1 is given in Table I. Note that they all differ only 
within the statistical error in the determination of the fractal dimension. We have verified that this result also applies 
to other different regions of parameter space, showing that this is not a coincidence for the instability bands shown. 
This is an important result since for every parameter region chosen, they will all have similar dynamical properties 
with each instability band then sharing similar properties. In this case, the dynamics is much less affected by the 
parameter space but only model dependent. In special, the fractality of those borders shows that one cannot be 
certain about the precise value of the Floquet exponent for a given mode, even though it is well inside an instability 
band. The numerical results are shown in Table |, for the first three bands and fo r thr ee different energy cutoff values. 



The results shown were produced considering = = 1 in Eqs. (2.29) and (2.37) 



We next study the infiuence of the number of field modes coupled to the infiaton field by changing the number of 



fermion fields and the boson fields ones N^. As it can be seen from Eqs. ( 2.3C ) and ( 2.38 ) this will, consequently, 
lead to changes to the dissipation coefficients. As we increase the decaying modes to both $ and Xj fields we expect 
the larger be the dissipation terms. The damping of the oscillations of the background field will eventually lead to 
a destruction of the instability bands. At the same time, even for underdamping oscillations we expect that as we 
increase the decaying modes the smoother will get the borders between stable and unstable regions of parameter space. 
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leading to a continuous decrease of the fractal dimension associated to the borders of the bands. This is confirmed by 
our numerical studies and the results for the change in fractal dimension as we increase the number of fields is shown 
in Fig. H, where we have restricted ourselves to the first instability band shown in Fig. 1, for an energy threshold of 
E = lOOi^o, and considered for simplicity = = N. Similar results would also follow to the other bands but are 
not shown here. 



IV. CONCLUSIONS 



Our results indicate that the fractal dimension of the borders (within the statistical error) does not change in any 
significant way for each instability band. It also shows very little dependence on the energy threshold chosen. These 
results, besides indicating that the border between the regions (defined by their own Floquet exponents) inside a given 
band is itself fractal and thus exhibits chaotic dynamics, shows that the dynamical properties of each instability band 
is independent of the precise region taken in the parameter space, and it is only model dependent. 

The fractal structure is also present if we fix the parameters k and g and vary only the initial inflaton conditions 
{if and ip). This shows a strong dependence of the outcome on the initial conditions, with major implications in 
the early universe models. For instance, an initial phase of warm inflation ||l6|, |2^, even if it is not able to generate 
enough radiation to match the standard Big Bang model, may turn out to be crucial for generating the adequate 
initial conditions for the prehe atin g phase. 

One may notice from Eq. (3J) that the effective mass squared m^fj = + fn^. + gj^{t) for the field Xk may 
become negative, triggering its decay via the so-called tachyonic instability (also known as spinodal decomposition). 
We have restricted the initial value of the inflaton so that it does not happen in our simulations. Nevertheless, the 
combination of such a decay with parametric amplification is an interesting and complex subject on its own and it 
will be addressed in a forthcoming publication |24]. 
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E 


4 


dV 


dV 


Eq 


3.26 ±0.06 


3.28 ±0.09 


3.22 ±0.07 


lOEo 


3.23 ±0.08 


3.23 ±0.09 


3.24 ±0.05 


100^0 


3.20 ±0.05 


3.20 ±0.07 


3.16 ±0.06 



TABLE I: Fractal dimensions df for the first (I), second (II) and third (III) instability bands, at different energy thresholds 
(E). Bands are numbered according to increasing , i.e., from left to right in Fig. |l| 
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FIG. 1: Band structure in the x space. Dark crosses indicate final energy above the cutoff (which, in this case, equals 
100 Bo). 
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FIG. 2: Fractal dimension d/ of the first band shown in Fig. 1 versus the number of interacting fields N . 



